Microcanonical temperature for a classical field: application to Bose-Einstein condensation 
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We show that the projected Gross-Pitaevskii equation (PGPE) can be mapped exactly onto Hamilton's equa- 
tions of motion for classical position and momentum variables. Making use of this mapping, we adapt tech- 
niques developed in statistical mechanics to calculate the temperature and chemical potential of a classical Bose 
field in the microcanonical ensemble. We apply the method to simulations of the PGPE, which can be used to 
represent the highly occupied modes of Bose condensed gases at finite temperature. The method is rigorous, 
valid beyond the realms of perturbation theory, and agrees with an earlier method of temperature measurement 
for the same system. Using this method we show that the critical temperature for condensation in a homo- 
geneous Bose gas on a lattice with a UV cutoff increases with the interaction strength. We discuss how to 
determine the temperature shift for the Bose gas in the continuum limit using this type of calculation, and obtain 
a result in agreement with more sophisticated Monte Carlo simulations. We also consider the behaviour of the 
specific heat. 



I. INTRODUCTION 

The Gross-Pitaevskii equation (GPE) has proven to be an 
extremely useful description of macroscopic Bose-Einstein 
condensates (BECs) at or near zero temperature (]]]. It is the 
first and sometimes only tool to be used in the description 
of many experiments in the field of non-linear atom optics 
and Bose-Einstein condensation. The validity of the GPE for 
many wide-ranging experimental situations now appears be- 
yond doubt. 

However, it has been proposed that the GPE can also be 
used to represent the non-equilibrium dynamics of Bose gases 
at finite temperature 01110 El- The underlying argument is 
that for modes of the gas with an average occupation much 
larger than one, the classical dynamics is far more important 
than the quantum dynamics. This is analogous to the semi- 
classical approximation utilised in laser physics for the elec- 
tromagnetic field. A major advantage of using the GPE in 
such situations is that it is non-perturbative, and so can be ap- 
plied in the region of the critical point as long as the condition 
on occupation numbers is observed. In Ref. |6j a finite tem- 
perature Gross-Pitaevskii equation is derived from the quan- 
tum many-body Hamiltonian for the Bose gas with this ap- 
proximation in mind. An alternative route to similar equations 
of motion is possible via the use of the Wigner representation 
0| . This approach may be more familiar to those from the 
quantum optics community. 

Some of the first numerical calculations utilising the GPE 
for finite temperature simulations were performed by Damle 
et al. js| and Marshall et al. @|. More recently there have 
been several calculations using the so-called "classical field" 
approximation. In particular we mention those of the Warsaw 
group the ENS group 
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and the 
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current authors lfl7i fisll . These calculations are all related by 
the fact that they include no damping terms in the GPE, and 
thus rely on ergodicity for the system to thermalize. Classical 
field methods involving both damping and stochastic terms 
have been considered by Gardiner et al. II 911. Stoof and Bi- 
jlsma HOIl , and Duine and Stoof PUl . 

While the qualitative results from previous work have been 
promising, there has been some difficulty in performing quan- 
titative calculations using such methods, and in particular in 
determining the temperature of the system at equilibrium. We 
partially addressed this issue in previous work using a variety 
of methods to determine the temperature of our simulations 
The most reliable of these involved fitting time- 
averaged quasiparticle occupations and energies to the clas- 
sical limit of the Bose-Einstein distribution function. How- 
ever, this method relies on the existence of a basis which ap- 
proximately diagonalizes the Hamiltonian (quasiparticles) for 
which energies and wave functions can be calculated in ad- 
vance. This method is therefore only applicable in the realm 
of perturbation theory, and fails for even moderate tempera- 
tures in systems with large nonlinearities. Hence it is desirable 
to find a more widely applicable scheme for unambiguously 
determining the temperature of numerical simulations. 

In a succinct yet insightful paper, Rugh l22ll expressed the 
temperature of a classical Hamiltonian system in terms of a 
phase space expectation value of a suitable function of the 
canonical position and momentum coordinates [Eq. dl9l of 
this paper.] Using the ergodic theorem, this expectation value 
over phase space can be interpreted as a dynamical average 
for a system in equilibrium, and immediately lends itself to 
application in numerical calculations. Rugh developed this 
procedure further in l23ll . and generalised it to include systems 
with other conserved quantities in addition to the energy l24ll . 
This generalization turns out to be crucial for the application 
of the method to the interacting Bose gas. Rugh's formula for 
the temperature has been applied to several systems to date, 
for example in the field of molecular dynamics. This has led 
to the notion of a configurational temperature for gases, which 
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only depends on the spatial co-ordinates of the particles, in 
addition to the usual kinetic temperature which only depends 
on the momenta lEll26ll27ll . 

In this paper we apply the microcanonical formalism of 
Rugh to the BEC Hamiltonian to determine the temperature 
of numerical simulations of thermal Bose gases. The method 
is non-perturbative and does not rely on the existence of well- 
defined Bogoliubov quasiparticles. The paper is organised as 
follows. In Sec.|n]we briefly summarise the expression for the 
temperature and other derivatives of the microcanonical en- 
tropy, and describe their application to the BEC Hamiltonian. 
Section ITTT1 presents our numerical results for our projected 
GPE (PGPE) system, while Sec. II VI relates our calculations 
to other dynamical calculations of classical (f> 4 field theory, as 
well as to calculations of the shift of the transition temperature 
for a homogeneous Bose gas. We conclude in Sec.lVl 



where can be chosen to be any scalar value, including zero. 
The vector field X can also be chosen freely within the con- 
straints 

VH ■ X = 1, VFi ■ X = 0, 1 < i < m. (5) 

Geometrically this means that the vector field X has a non- 
zero component transverse to the H = E energy surface, and 
is parallel to the surfaces Fj = Ii. The expectation value in 
Eq. Q is over all possible states in the microcanonical en- 
semble; however if the ergodic theorem is applicable then it 
can equally well be interpreted as a time-average. For fur- 
ther details on the origin of this expression we refer the reader 
to Rugh's original papers I2 3 . [23ll24ll . as well as derivations 
found in Giardina and Levi l28llT .Tepps et al. l26ll and Rick- 
ayzen and Powles l29ll . 



II. FORMALISM 



B. Dimensionless BEC Hamiltonian 



A. Hamiltonian 

We consider a classical system with M independent modes. 
The Hamiltonian can be written as H = H(T), where T = 
{Ti} = {Qi, Pi} is a vector of length 2M consisting of the 
canonical position and momentum co-ordinates. In these co- 
ordinates we define the gradient operator V in terms of its 
components V, = d/dTi. 

In the notation of Rugh ll24ll . the Hamiltonian H may have 
a number of independent first integrals, labelled by F — 
Fx, . . . , F m , that are invariant under the dynamics of H. We 
could define Fq = E, and include the conserved energy with 
the other constants of motion in this notation, but for clarity 
we consider it separately. A particular macro-state of such a 
system can be specified by the values of the conserved quan- 
tities, labelled as H = E,Fi = Ii. 

The expression for the temperature of such a system in the 
microcanonical ensemble is given by 
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(1) 



where all other constants of motion are held fixed, and where 
the entropy is given by 

e s/kB = JdT S[E-H(T)] Jp[/; - F^T)]. (2) 



In this case, the temperature of the system can be written as 
1 



i>-x(r) , 



(3) 



where the angle brackets correspond to an ensemble average, 
and the components of the vector operator T> are 
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The full quantum many-body Hamiltonian for the Bose gas 
in dimensionless form is 

H = J d 3 x V* t (i)-V*(i) + V r (i)* t (i)*(i) + 

-* t (x)# t (x)*(x)*(x) , (6) 



where H = Ne^H, N is the number of particles in the sys- 
tem, x = x/i, L is the unit of length, = h 2 / (2mL 2 ) is the 
unit of energy, m is the mass of the particles and V(x) is the 
dimensionless external potential if any is present. The dimen- 
sionless quantum Bose field operator ^(x) is here normalized 
to one, J d 3 x(>f ^(x)'I'(x)) = 1, and C n \ is the nonlinear con- 
stant defined as 



C n l 
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e L L 3 



(7) 



where a is the s-wave scattering length. In this expression 
we have assumed a high momentum cutoff and made use of 
the replacement of the true interatomic potential with the two 
body T-matrix, V(x) — > Uq5(-k), where Uo = <lTrh 2 a/m. 

In Ref. 01 the field operator is split into a classical part 
and a quantum part, with the boundary determined by the re- 
quirement that the average occupation number (Nk) of modes 
below the cutoff satisfies (Nk) ^> 1. Equations of motion 
were derived for the classical part, before taking the mean 
value. This resulted in the finite temperature Gross-Pitaevskii 
equation (FTGPE), which describes the evolution of a classi- 
cal field coupled to an effective bath described by a quantum 
Boltzmann-like equation. This equation proves to be some- 
what difficult to solve numerically, and in Refs. fl7lfl^l we re- 
ported results focussing on a simplification we termed the pro- 
jected Gross-Pitaevskii equation (PGPE). This equation de- 
scribes the evolution of a classical field only, with a cutoff at 
a given momenta or energy. It is identical to the usual GPE 
except that it evolves a wave function which is restricted to a 
finite-sized basis satisfying the classical condition (Nk) 3> 1. 
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The PGPE for a homogeneous system is written explicitly for 
the homogenous gas in Eq. J26I of this paper. 

In this paper we wish to determine the temperature of the re- 
stricted system described by the evolution of the PGPE. Thus 
the Hamiltonian we consider is the classical version of Eq. (|6} 
obtained by replacing the field operator ^(x) with the classi- 
cal field ^(x), subject to the important restriction that ^>(x) 
is constructed from a finite number of low-energy modes. We 
can therefore write it in the form 



(8) 



where C labels the classical modes in the coherent region be- 
low the cutoff, as defined in @|. 



C. Canonically conjugate position and momentum variables 

We must now make a choice of the canonically conjugate 
co-ordinates of our Hamiltonian. As we are defining our clas- 
sical field in a basis, it seems natural to convert to a basis 
representation. If we choose our basis to be that which di- 
agonalises the ideal gas Hamiltonian [the first two terms of 
Eq. l|6}] we find 



H 



/ j e nC n C n 



iClc p c q (mn\pq), (9) 



ninpq 



where the matrix element is 



(mn\pq)= / d 3 x^(x)0;(x)</> p (x)^ g (x). 



(10) 



The equation of motion for the {c n } is given by the PGPE. 
This problem can be mapped exactly to the one considered by 
Rugh by defining real, canonically-conjugate coordinates Q n 
and P„ 



^=« + c„), Pn=i\h!r{c* n -c n ), (11) 



with the corresponding inverse transformation 





<2en 

(12) 

With these definitions, the evolution of the c„ coefficients 
given by the PGPE maps exactly to the evolution of the co- 
ordinates Q n and P n given by Hamilton's equations. The 
PGPE is therefore in one-to-one correspondance with a classi- 
cal microcanonical system and its equilibrium propreties can 
be studied using the wide variety of techniques which have 
been developed in classical statistical mechanics. 

We have performed numerical calculations for the homoge- 
neous PGPE, and so we use a plane wave basis where c/> n (x.) 
= exp(zk„ -x), n = {n x , n y , n z } and e„ = \k n \ 2 = (27r|n|) 2 . 
However, the method we describe is general and can be ap- 
plied directly to inhomogeneous systems for BECs in mag- 
netic and optical dipole traps. 



Calculations on a grid 

The implementation of a projection operator in the GPE is 
an essential feature of any classical simulation. While we have 
explicitly defined a projection operator in terms of a basis set, 
other authors have implicitly chosen a momentum cuto ff by 
the use of a finite-size grid in their GPE simulations fulfill . 
The method of temperature determination described in this pa- 
per can also be applied to these calculations, but with a differ- 
ent choice of postion and momentum co-ordinates. 

On a finite grid, the Hamiltonian © can be discretised in 
real space and the classical equivalent can be written as 



H = h x h y h z ^ 



(Va n ) 2 + (V/3„) 2 + 



(13) 



5 } labels the grid point, h x , h y , h z are 



where n = { 

the grid spacings for each axis, and we have defined ip n = 
a n + ij3 n . In this case the appropriate position and momentum 
variables are 



Q n = V2a n , P n = V20 n 



(14) 



However, we believe that it is important to define the pro- 
jector using a basis that is relatively well-defined in energy at 
the cutoff (we stress that this does not mean that the basis has 
to be well-defined in energy below the cutoff). It has previ- 
ously been shown ll30ll that the single particle energy levels of 
a partially condensed system are essentially those of the trap- 
ping potential for energies e > En « 3/ics where [ic is the 
condensate energy eigenvalue. Thus the above cutoff projec- 
tor can be written 



Q{^(x)} = 



0fc(x) 



:(x')F(x'), (15) 



where the {<frk} are the basis states appropriate to the poten- 
tial, and the notation k ^ C describes a summation over all 
modes above the energy cutoff Er. As this basis is complete, 
the below cutoff projector is simply 



v = i-q, 



(16) 



which gives the result written explicitly in Eq. d27> . We also 
require the classical condition 



k B T 
E R - \i 



» 1, 



(17) 



to hold at the cutoff and so for Er w this should also be 
satisfied. 

For a trapped Bose gas, the implicit momentum projector 
based on the finite-grid method is not at all well-defined in 
energy at the cutoff, and we believe that this may lead to diffi- 
culties. However, this is yet to be invesigated numerically; for 
further discussion of this issue we refer the reader to 13111 . 
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D. Choice of vector field X 

In order to satisfy the conditions we can choose a vector 
field of the form 



X = aVH + ^2biDFi, 



(18) 



where the m + 1 coefficients {a, hi} are determined by the 
m + 1 simultaneous equations in Eq. 0. Due to the freedom 
in the choice of the vector operator V we can set any compo- 
nent of the length 2M vector X to zero. This turns out to be 
useful as the components corresponding to the momentum and 
position variables can be different orders of magnitude. Two 
particular choices we make use of later are Xp with T> = 
V P = {0,8/dPi} and X Q with V = V Q = {d/dQ ll 0}. 
These lead to two different calculations for the temperature 
that only agree in general if the system is in thermal equilib- 
rium. This provides a useful check that the simulations have 
in fact thermalized. 

In Rugh's first two papers ll22l l23ll the only first integral 
considered was the energy, and he chose T> = V which 
yielded the (dimensionless) formula 



1 

T 



V 



VH 

WW 



(19) 



itly in 12411 for a system of particles with a conserved centre- 
of-mass motion. 

An exception to this, however, is the conservation of nor- 
malization N = J2 n c n c n- This must be considered explicitly 
because there is no co-ordinate system in which it can be made 
to vanish. The constraint on TV means that the ground state of 
the system will, in general, have a finite energy. For exam- 
ple, a non-interacting gas in a harmonic trap of frequency uj 
must have at least the zero-point energy hoj/2 for each spa- 
tial degree of freedom. For a non-ideal, homogeneous gas the 
restriction that at least one of the c n must be non-zero means 
that there will always be a finite interaction energy associated 
with the ground state energy Eq = C n \/2. These energy con- 
tributions are not accessible for thermalization, however, and 
including the normalization constraint allows them to be re- 
moved. We note, however, that the effect of this constraint is 
in general more complicated than a simple subtraction of the 
ground state energy (which could be achieved by hand) and 
depends on the definition of the operator V used to calculate 
the temperature, as shown below. 

To deal with the normalization constraint, we need to 
choose a vector field X which satisfies Eqs. © with F\ = 



n C n C n 



The result is 



X = 



VH - \ N VN 
VH\ 2 - \ N VN -VH' 



(20) 



For the BEC Hamiltonian we consider, however, there are 
other first integrals that must be taken into account. Most 
importantly, the evolution conserves the normalization of the 
wave function, but other first integrals that may occur are both 
the angular and linear momentum. 

The effect of including these additional first integrals in the 
definition of the vector field X is to account for the energy that 
is associated with a conserved quantity and hence is unavail- 
able for thermalization. This ensures that only the appropriate 
free energy is used to calculate the temperature. We conjec- 
ture that the same result can be achieved by first transforming 
to a co-ordinate system where the total angular and linear mo- 
menta, etc, are all zero and therefore do not contribute to the 
energy of the system. In fact, Rugh demonstrated this explic- 



where the parameter 



VN -VH 
\VN\ 2 ' 



(21) 



looks similar to a chemical potential. For a system on a real 
space grid with T> = V and a Hamiltonian given by Eq. Jl 3I > 
we find that \m = Mope, where /Ugpe is the usual Gross- 
Pitaevskii form of the chemical potential, obtained from the 
Hamiltonian of Eq. d 1 3 b by doubling the interaction term. 
However, in general the expression of Eq. Mil does not have 
a simple interpretation. 

Substituting Eq. d20b into Eq. Q we find that our full ex- 
pression for the temperature is 



J 



1 

T 



V 2 H - X N V 2 N - VX N ■ VN 
\VH\ 2 -\ N {VH -VN) 



UVH - \ N VN) ■ [V\VH\ 2 



(VH ■ VN)V\ 



N 



VH\ 2 - \ N (VH -VN)f 



X N V{VH-VN)] ^ _ (22) 



The second term in this expression is of order 1/M, and so in 
many situations it can reasonably be neglected. However we 
have calculated the full expression for all results presented in 
this paper. 



E. Other thermodynamic quantities 



The method described in this paper can also be used to cal- 
culate first derivatives of the microcanonical entropy with re- 
spect to any of the first integrals of the Hamiltonian l24ll . In 
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particular, we find to calculate the quantity 



ur 3 / E.Fi 



the constraints on our vector field should be 

VH-X = 0, VF l -X = l, VFj-X = 0, i^j. (23) 
For the BEC Hamiltonian, we have 



dS 
~dN 



(24) 



and implementing the required constraints, we find an appro- 
priate vector field is that given by Eq. (12012 1> but with the 
roles of H and N reversed. 

In addition, higher order derivatives of the entropy can also 
be determined, making available quantities such as the spe- 
cific heat c sp of the system l23ll . This quantity could in prin- 
ciple be calculated from the expression 



1 



= 1 - 



• S|> 



(V ■ (XV ■ X)) 
(V • X) 2 



(25) 



where the vector X is determined by Eq. 112012 It . However, 
for the BEC Hamiltonian the expressions for such quantities 
are unreasonably complicated, and we do not consider them in 
this paper. Instead, higher derivatives will simply be obtained 
numerically once the temperature is determined. 



in. NUMERICAL RESULTS 



In this section we apply the new formula Eq. \22 \ to data 
from simulations of the PGPE described in Il7l Hal , as well 
as to many new simulations with a wider range of energies 
and nonlinear parameters C n \. For a full description we refer 
the reader to Ref. fiaV Briefly, the calculations evolve the 
projected Gross-Pitaevskii equation [;6] for the homogeneous 
gas in three dimensions 



% dr 



-V^(x) + C nl P{|7/.(5c)|>(5c)}. 



(26) 



The nonlinear constant is C n i = 2mNUo/h L, where N is 
the total number of particles in the volume, and L is the period 
of the system. Our dimensionless parameters are x = x/L, 
wave vector k = kL, energy e = s/sl, and time r = ejj.j'h, 
with el = ft 2 /(2mL 2 ). The projection operator V excludes 
all components of the nonlinear term in the GPE outside the 
coherent region, and is defined by [c.f. Eqs. ( I15I16H 



V{F(x)} = 



kec 



d 3 x' ^.(x')F(x'), (27) 



where {<j>k} is an orthonormal basis appropriate to the prob- 
lem. For the homogeneous system with periodic boundary 
conditions the relevant basis is the plane wave states, and so 
this procedure is simply the application of a forward fourier 



transform, removal of components with k > k c , followed by 
the inverse transformation. The quantity k c defines the mo- 
mentum cutoff for the coherent region, and for all data pre- 
sented in this paper we use k c = 15 x 2tt. 

We begin with randomised initial fields ip(x) with a given 
energy on a 3D grid with 32 points in each dimension, and 
evolve these until the field has reached equilibrium. We calcu- 
late all thermodynamic quantities from sampling two hundred 
field configurations in equilibrium. 



Cutoff dependence of simulations 

The choice of momentum cutoff k c = 15 x 2ir is motivated 
simply by computational convenience. It also allows for com- 
parison of the Rugh method of temperature measurement with 
data from earlier calculations. 

For a given initial energy, the resulting equilibrium tem- 
perature depends on the number of modes below the cutoff. 
This can be easily understood from the equipartition theorem 
— if more modes are present, less energy will be contained 
in each one and therefore the final temperature will be lower. 
Also, the dimensionless critical temperature for a system with 
a fixed normalisation depends on the cutoff, as can be seen in 
the text beneath Eq. J28I . 

Work is currently in progress to develop a description of the 
modes above the cutoff and their coupling to the PGPE. The 
aim of this work is a complete computational method which 
will be insensitive to the exact position of the cutoff. Explor- 
ing and developing techniques for the non-perturbative classi- 
cal field is an important part of this programme, and we focus 
on this aspect of the problem in this paper. 

Despite this, there are some equilibrium calculations which 
can be carried out immediately using an approximate treat- 
ment of the modes above the cutoff. We present results for 
one such calculation (the shift in T c with interaction strength) 
in this paper. These results have only a weak dependence on 
the cutoff. 



Use of the classical field method at T c 

The classical field can only describe modes that satisfy the 
high occupation condition. But even at the critical tempera- 
ture and above, the lowest energy states will have the largest 
occupations — and for a wide range of parameters, many of 
these can satisfy Nk S> 1. These are the modes that are re- 
sponsible for critical behaviour, such as the shift in T c and the 
increase in specific heat. The remaining modes (that are not 
simulated) behave essentially as an ideal gas. 

As a physical example, consider our simulations for C n \ = 
20000. Choosing L = 25/j.m, and Rb-87, this corresponds 
to approximately 3.8 x 10 5 atoms below the cutoff satisfying 
N k > 10 at T c of about 370 nK. There are about 1.3 x 10 6 
atoms in total, with a total number density of 8.3 x 10 13 cm -3 . 
Thus in this situation nearly 30% of the atoms are simulated 
by the PGPE. 
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FIG. 1 : Plot of the relative differences of simulation temperatures Tq 
and Tp calculated from Eq. i22\ with temperatures To determined 
from the same data plotted in Fig. 9 of Ref. Hill . Grey triangles: 
AT = T P - T , black dots: AT = T Q - T . (a) C n i = 500, (b) C n i 
= 2000, (c) C* n i = 10000. 



A. Comparison of methods of temperature determination 

As described in Sec. lllDl there are many choices of the op- 
erator T> that may be used in Eq. \22\ . The resulting calcula- 
tions only give the same temperature if the system is in equi- 
librium, so this provides an important confirmation that the 
system has thermalized. In this paper we consider two cases 
V Q = {d/dQi, 0} and V P = {0, d/dP l }, and we refer to 
the temperatures calculated from these operators as Tq and 
Tp respectively. Allowing Q or P derivatives only in the sep- 
arate definitions of the operator T> simplifies the calculation of 
Eq. (I22> due to the elimination of mixed derivatives. 

We begin by comparing Tq and Tp with previous results 
from Ref. IU8I1 . In this earlier work we obtained temper- 
atures using three different methods, two based on Bogoli- 
ubov quasiparticles and perturbation theory and a third non- 



perturbative calculation. This third method did not have a 
firm theoretical basis; however, we showed that the results 
were consistent with the two other calculations in their regime 
of validity, and gave reasonable results more generally. Fig- 
ure [2 shows the relative differences between the new simu- 
lation temperatures Tq and Tp calculated from Eq. i22\ and 
the temperatures To determined from the earlier method three. 
The simulation data used is the same as that plotted in Fig. 9 
of Ref. O. 

We can see from Fig. [2 that only a small number of points 
differ by more than one percent from the previously deter- 
mined values, and even these would be hard to detect on a 
plot of the absolute temperatures. These results therefore val- 
idate our earlier non-perturbative method for temperature de- 
termination in a homogeneous system (described in Sec. VI.D 
of Ref. 11811 '). Figure ^ also shows that in general the values 
of Tq and Tp agree with each other within their error bars. 
The error is determined from the standard deviation of the ex- 
pectation value of Eq. ( 1221 divided by the square root of the 
number of samples (in this case two hundred). This estimate 
assumes gaussian statistics, which seems reasonable when the 
distribution of values is plotted as a histogram; however, it 
may underestimate the actual error somewhat. The agreement 
between these distinct determinations of temperature confirms 
their validity and provides important further evidence that the 
PGPE evolves randomised initial states to a thermodynamic 
equilibrium consistent with the microcanonical ensemble. 



B. Shift of the transition temperature 

Figure|3a) plots the equilibrium temperatures and conden- 
sate fractions for several series of simulations with different 
nonlinearities C n \, as well as for the ideal gas. These can be 
interpreted to be simulations at a fixed density with a varying 
scattering length. It is immediately obvious that qualitatively 
the transition temperature increases with increasing nonlin- 
earity, and this was noted in fl8ll . Much more data has been 
collected for this paper, and we now have a much more reli- 
able measure of temperature. Thus we can now look at the 
shift of the critical temperature with the nonlinear parameter 
C n i for our PGPE system. 

We can calculate the transition temperature for a non- 
interacting gas with equipartion occupation numbers and a 
momentum cutoff k c in the continuum limit via 



N = N + 



k < d 3 k 



(2tt) 3 h 2 k 2 /(2m)-n' 



(28) 



We find that the dimensionless critical temperature for a ho- 
mogeneous PGPE system with a momentum cutoff of k c = 
2-kk/L is T c (C n i = 0) = tt/k, where the dimensionless tem- 
perature is defined by T = ksT / {Nej). 

Identifying the critical point in a finite-sized system with 
interactions, however, is somewhat more difficult. Here we 
make use of the method of Binder cumulants 1 32], which have 
been used in other finite-size calculations for the Bose gas 
l33ll . We note that the theory behind Binder cumulants is de- 
rived entirely from canonical statistical mechanics. However, 
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curves of C vs T intersect for different lattice sizes is univer- 
sal for a given universality class, which is three-dimensional 
XY for our system. It has been calculated by Campostrini et 
dl. Q that this critical value is C c = 1.2430(1)(5), where 
the first number in parentheses is due to statistical errors and 
the second is due to systematic errors. 

We therefore determine the critical temperature from our 
simulations by finding the energy at which the Binder cumu- 
lant takes the value C c in equilibrium. Due to our limited 
statistics from 200 field samples, the results are somewhat 
noisy, but we are able to identify T c for the simulations to 
an accuracy of approximately one percent. 

We note that for the case of C n \ = 20000, the predicted 
shift in critical temperature is more than 60%. However, this 
corresponds to the shift in dimensionless temperature of the 
low-energy states, not the shift in the critical temperature of 
the complete system which will be smaller. This will be dis- 
cussed in more detail in Sec. II V Bl 




FIG. 2: (a) Plot of the condensate fraction versus temperature for a 
number of interaction strengths. Solid line C n \ = 0, crosses C n i = 
500, solid dots C n i = 2000, open circles C n i = 5000, plusses 
Cni = 10000, stars C n i = 15000, open triangles C n i = 20000. 
(b) Plot of the transition temperature versus interaction strength. The 
transition temperature is determined by the method of Binder cumu- 
lants as described in the text. 



the calculations of Caiani et al. I34i 13511 suggest that it is valid 
as a numerical tool in the microcanonical ensemble, and we 
shall follow their lead. The Binder cumulant can be written as 



C = 



(No)' 



(29) 



where Nq is the population of the zero momentum conden- 
sate mode in our simulations. This quantity changes smoothly 
from one for the condensed system (ordered phase) to two for 
the uncondensed system (disordered phase), with the width 
of the transition region decreasing with increasing lattice size. 
However, in lattice field theory the chemical potential at which 



C. Calculation of the specific heat 

Although the specific heat can theoretically be determined 
by a similar procedure to that used for the temperature, the ac- 
tual formulae are rather complicated and difficult to calculate. 
Instead in this section we use numerical methods to calculate 
curves for the specific heat. 

The calculation of numerical derivatives is difficult for data 
with statistical errors. Here we have used a smoothing spline 
fitting technique to the raw numerical data for energy and tem- 
perature, and calculated the derivative from this fit. Examples 
of the spline fits to the numerical data are plotted in Fig. [3] 

The specific heat curves calculated from the data in Fig. [3] 
are shown in Fig.Ua). The units of the vertical axis are scaled 
by the specific heat of the ideal Bose gas for the same system 
at T = 0. We can see that there is a strong peak near the crit- 
ical temperature that increases with increasing C n \. Scaling 
theory for critical points in the thermodynamic limit suggests 
that the specific heat will be discontinuous at the phase tran- 
sition. In our case the peak is not exactly at T c , as perhaps 
would be expected. We presume that this is due to to a combi- 
nation of finite size effects and numerical errors in the fitting 
procedure, which we estimate to be a few percent. Similar 
behaviour has also been noted in l34ll . Figure|4jb) shows the 
maximum value of the specific heat plotted versus C n \. 



IV. RELATION TO OTHER WORK 

A. Dynamics of </> 4 lattice field theory 

The results presented in this paper for the homogeneous 
Bose gas have many similarities to classical </> 4 lattice field 
theory, which is often studied in relation to second order phase 
transitions. In such studies the field is discretised on a lattice 
with the spatial derivatives of the Hamiltonian being approxi- 
mated by finite difference methods. Monte Carlo simulations 
are then performed to study the thermodynamics. 
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FIG. 3: The system energy plotted against the temperature, scaled in 
units of the critical temperature T c so that the curves are distinct in 
the linear region. Solid line C n \ — 0, solid dots C n \ — 2000, plusses 
C n i = 10000, open triangles C n \ = 20000. The grey lines are the 
smoothing spline fits to the numerical data. 



However, there have also been "molecular dynamic" sim- 
ulations of such field theories, and in particular we note the 
work of Caiani et ai, who have considered the phase tran- 
sition via dynamical simulation of the </> 4 model in both two 
ll35ll and three dimensions ll34ll . Their equations of motion 
are distinct from those of this paper by virtue of being second 
order in time. Their paradigm Hamiltonian in <i-dimensions is 



#[</>] = 



1 



J 



X — 7T 

2 



(x) + -[V#x)]' 



1 



A 



2 -0 2 (x) + I 4 (x), 
(30) 

with the canonical position variables </>(x) and conjugate mo- 
menta 7r(x) = 0(x), where is a vector quantity with up 
to four dimensions. We note that as this Hamiltonian is of 
the form H = 7r 2 /2 + V((j>), both the temperature and spe- 
cific heat of these simulations can be calculated from expec- 
tation values of the kinetic energy. This is not possible for 
the Hamiltonian we consider in this paper where the interac- 
tion term mixes powers of the position and momentum co- 
ordinates. 

Also, in Ref. 13511 the parameters used were J = 1, A = 
0.6, and for Ref. yj the values J = 1, A = 0.1 and A = 4 
are specifically mentioned. Thus these calculations appear to 
be in quite a different regime to the results presented here. 
Despite these differences, however, it seems that much of their 
numerical data is qualitatively similar to ours. 



B. Shift of T c in the continuum limit 

The results presented in this work can also be connected to 
the issue of the shift in the transition temperature for the ho- 
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FIG. 4: (a) The numerically calculated specific heat curves for var- 
ious interaction strengths. The peaks occur at temperatures a few 
percent below the identified transition temperature. We estimate the 
error for these curves to be of order a few percent. Solid black line 
Cni = 0, dashed line C n i = 2000, dotted line C n i = 10000, solid 
grey line C n \ = 20000. (b) The maximum value of the specific heat 
plotted versus the dimensionless interaction strength C n i. For both 
(a) and (b), the specific heat is plotted relative to the corresponding 
value at T = 0, and so the quantities are dimensionless. 



mogeneous Bose gas, which has been the subject of a number 
of recent papers. In the weak interaction limit the shift AT C 
has the form 



Tco 



= can 1 '*. 



(31) 



where T c q is the transition temperature for the ideal gas, n 
is the density, a is the s-wave scattering length and c is a di- 
mensionless constant. The value of c cannot be determined by 
perturbation theory as this breaks down at second-order phase 
transitions due to infrared divergences. There have been sev- 
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eral calculations of the value of c, differing by up to an order 
of magnitude and even in sign (see the summary in l37io . 

The dimensionless constant c has recently been determined 
via Monte Carlo calculations by Arnold and Moore l33l Ijl 
and Kashurnikov et al. 13811 to be c = 1.32 ± 0.02 and 
c = 1.29 ± 0.05 respectively. These calculations were carried 
out via classical 4 field theory, which can be systematically 
matched to the problem of the homogenous interacting Bose 
gas. The Monte Carlo calculations proceeded by sampling the 
classical action 



S 
1 



d 3 > 



V*(x) - 



(32) 

on a lattice at a fixed temperature T, where (3 = (fc^T) -1 . 
The value of /j, c g was adjusted until the critical point was 
reached, and thus the shift in critical density n c = (\tp\ 2 ) from 
the ideal gas value n C Q could be measured. The shift in critical 
temperature at a fixed density can then be determined from 



2 An r 



T, 



(•() 



3 n c0 



(33) 



which is easily derived from the formula for the critical tem- 
perature of the ideal gas. While this procedure seems straight- 
forward, in practice it is necessary to give careful considera- 
tion to finite-size effects in the calculation — see Ref. l33ll for 
a detailed discussion of these matters. 

The results of simulations similar to those presented here 
can also be used to calculate a value for c, as we are also 
sampling the thermodynamic functions of classical </> 4 field 
theory. The Monte Carlo calculations fix the temperature and 
adjust the value of fi c g in Eq. d32i which then determines the 
normalisation of the field. In our calculations, we adjust the 
energy of the initial state to find the critical point and deter- 
mine the temperature using the method described above. Our 
simulations have a fixed normalisation, but the dimensionless 
temperature T oc T/N, so for a given value of C n \ we can in- 
terpret our results as being at a fixed temperature and a varying 
density. 

The main difference between the methods is the manner 
in which field configurations are sampled. The Monte Carlo 
methods can use the most efficient update possible, as long as 
the samples are canonical at a given temperature. Our calcula- 
tions solve for the evolution of a microcanonical field, and use 
the theorem of ergodicity to generate an ensemble. We have 
one minor advantage in that our momentum cutoff is spheri- 
cally symmetric, whereas the Monte Carlo calculations simu- 
late the first Brillouin zone of the lattice. However, the molec- 
ular dynamics method suffers from critical slowing down — as 
the energy of the highest modes is oc fc 2 , we require time steps 
of order St = 1/fc 2 where k c is the momentum cutoff. Thus 
our simulations are disproportionally less efficient for larger 
grids compared to the Monte Carlo calculations, and will not 
be able to gene rate results as accurately for a given compu- 
tation time I39B . Nonetheless, we can use our simulations to 
confirm qualitatively the results of the Monte Carlo analysis, 
providing an independent demonstration of the validity and 
potential usefulness of our temperature determination. 




10 15 20 25 



10 4 an 1/3 



FIG. 5: Shift in the critical temperature with interaction strength de- 
termined from the results presented in this paper with A r beiow(C n i = 
0) = 10 10 . The dashed line is a linear fit to the first two data points 
and this has a slope of 1.3 ± 0.4. 



As a simplified illustration, we follow through the logical 
procedure that would be required to calculate a value for c. To 
consider the shift in the critical point, we can consider the shift 
in the critical density given a fixed critical temperature Tq. In 
our numerical simulations we keep C n \ = %iraNY, c \ ov , / L fixed 
and measure a shifted critical temperature 



^VbolowCL 



(34) 



Here iVbeiow is the number of particles below the cutoff. If 
we fix the critical temperature at To as well as the system size 
L (and hence el), we can interpret the increase in the dimen- 
sionless quantity T c as a decrease in the value of -/Vbeiow and 
hence a decrease in the critical density. The most important 
point to note is that as long as we have k c 3> fco, where fco 
labels the division between quasi-particle and particle like ex- 
citations at the transition temperature, then particles above the 
cutoff will not be significantly affected by the change in the 
interaction strength. 

We therefore calculate -/V a bovc = N to t — N^ c \ ow for the 
ideal gas, where iV to t is the total number of particles. This 
will be a constant as long as k c ^> fco. We can then calculate 
M>eiow (C n i) and hence the shift in the critical density from 
the simulation data, and by using the relation of Eq. i33\ . we 
obtain the shift in the critical temperature. 

This can then be plotted against an 1 / 3 and the slope at 
the origin determines the coefficient c. This plot is given in 
Fig. E] where we have set N hc i ow (C nl = 0) = 10 10 . We 
note that the method does not depend on the value chosen for 
^Vbdow(Cni = 0) as long as it is large enough that (Nk) ^> 1 
is well satisfied. 

By fitting a straight line to the first two points as illustrated 
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in Fig. |3 we get an estimate for the coefficient 

c= 1.3 ±0.4, (35) 

where the error specified is due to the uncertainty in the value 
of T c for the data poi nt. This agrees with the value deter- 
mined in Refs. 1371 13811 — a result that should be treated with 
caution. The correct value of c will only be reached in the 
limit of large volume and small lattice spacing, and we be- 
lieve we have not reached this regime. For comparison with 
the results of Arnold and Moore, for our first data point we 
have Lu w 325 and uai att ps 10.2, where u = 3TC n \/L and 
aiatt = L/32. Our other data points have values for these 
quantities that are much larger than this. Arnold and Moore 
suggest that Lu > 400 and ua\ at t < 6 are necessary to get 
an accurate result for c without a finite-size scaling analysis 

ElEi. 

We could potentially improve our results by performing 
such a finite-sized scaling analysis, but there is little reason 
to do so given the greater accuracy obtained in Refs. 137113811 . 
The purpose of this calculation is to demonstrate a useful ap- 
plication of our temperature determination with the PGPE in 
a non-perturbative regime. In this regard the qualitative agree- 
ment with earlier more involved and specialized calculations 
provides a pleasing confirmation of the general validity of the 
method. 



V. CONCLUSIONS 

We have shown that the projected Gross-Pitaevskii equa- 
tion can be exactly mapped to Hamilton's equations of motion 
for canonically conjugate position and momentum variables. 



Using this mapping we have described how to utilise the mi- 
crocanonical thermodynamic method of Rugh 12411 to measure 
the temperature of PGPE simulations in the non-perturbative 
regime. This method agrees with previous calculations de- 
scribed in Ref. fiUl . but has a rigorous theoretical justifica- 
tion and wider applicability. Using this approach, we have 
quantitatively measured the shift in the critical temperature 
for condensation with the nonlinear constant C n \. We have 
also observed that the specific heat reaches a maximum near 
the transition point as expected from the theory of continu- 
ous phase transitions, and that the peak value increases with 
the nonlinearity. Finally, we have made a connection between 
these calculations and Monte Carlo simulations that have de- 
termined the shift in the critical temperature with scattering 
length of the homogeneous Bose gas in the continuum limit. 
This is further evidence that the projected GPE should be valid 
for dynamical calculations through the critical region as long 
as the condition on the occupation numbers is satisfied. 



Acknowledgments 

The authors are grateful to Peter Arnold and Guy Moore for 
their insightful comments on the relation of this work to that 
of Refs. B33i 13711 . We also thank Crispin Gardiner for his use- 
ful comments. MJD would like to thank Tim Vaughan, Karen 
Kheruntsyan, Joel Corney, and Peter Drummond for several 
useful discussions at various stages of this work. MJD ac- 
knowledges the financial support of the University of Queens- 
land and the Australian Research Council Centre of Excel- 
lence grant CE0348178. SM would like to thank the Royal 
Society of London for financial support. 



[1] F. Dalfovo, S. Giorgini, L. P. Pitaevskii and S. Stringari, Rev. 

Mod. Phys. 71, 463 (1999), and references within. 
[2] B. V. Svistunov, J. Moscow Phys. Soc. 1, 373 (1991). 
[3] Yu. Kagan, B. V. Svistunov, and G. V. Shlyapnikov, Zh. Eksp. 

Teor. Fiz. 101, 528 (1992) [Sov. Phys. JETP 75, 387 (1992)]. 
[4] Yu. Kagan and B. V. Svistunov, Zh. Eksp. Teor. Fiz. 105, 353 

(1994) [Sov. Phys. JETP 78, 187 (1994)]. 
[5] Yu. Kagan and B. V. Svistunov, Phys. Rev. Lett. 79, 3331 

(1997). 

[6] M. J. Davis, R. J. Ballagh, and K. Burnett, J. Phys. B 34, 4487 
(2001). 

[7] M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, S. M. 
Tan, M. J. Collett, D. F. Walls and R. Graham, Phys. Rev. A 58, 
4824 (1998). 

[8] K. Damle, S. N. Majumdar, and S. Sachdev, Phys. Rev. A 54, 
5037 (1996). 

[9] R. J. Marshall, G. H. C. New, K. Burnett, and S. Choi, Phys. 
Rev. A 59, 2085 (1999). 
[10] K. Goral, M. Gajda, and K. Rzazewski, Opt. Express 8, 92 
(2001). 

[11] K. Goral, M. Gajda, and K. Rzazewski, Phys. Rev. A 66, 
051602 (2002). 

[12] H. Schmidt, K. Goral, F. Floegel, M. Gajda, and K. Rzazewski, 
J. Opt. B, 5, 96 (2003). 



[13] A. Sinatra, Y. Castin and C. Lobo, J. Mod. Opt. 47, 2629 (2000). 
[14] A. Sinatra, C. Lobo, and Y. Castin, Phys. Rev. Lett . 87, 210404 
(2001). 

[15] A. Sinatra, C. Lobo, and Y. Castin, J. Phys. B 35, 3599 (2002). 
[16] C. Lobo, A. Sinatra, and Y. Castin, cond-mat/0301628 
[17] M. J. Davis, S. A. Morgan, and K. Burnett, Phys. Rev. Lett. 87 
160402 (2001). 

[18] M. J. Davis, S. A. Morgan, and K. Burnett, Phys. Rev. A 66 
053618 (2002). 

[19] C. W. Gardiner, J. R. Anglin, and T. I. A. Fudge, J. Phys. B 35, 
1555 (2002). 

[20] H. T. C. Stoof and M. J. Bijlsma, J. Low. Temp. Phys 124, 431 
(2001). 

[21] R. A. Duine and H. T. C. Stoof, Phys. Rev. A, 65, 013603 
(2001). 

[22] H. H. Rugh, Phys. Rev. Lett. 78, 772 (1997). 

[23] H. H. Rugh, J. Phys. A, 31 7761 (1998). 

[24] H. H. Rugh, Phys. Rev. E 64, 055101 (2001). 

[25] B. D. Butler, G. Ayton, O. G. Jepps, and D. J. Evans, J. Chem. 

Phys. 109, 6519, 1998. 
[26] O. G. Jepps, G. Ayton, and D. J. Evans, Phys. Rev. E 62, 4757 

(2000) . 

[27] J. Delhommelle and D. J. Evans, J. Chem. Phys. 114, 6229, 

(2001) ; ibid. 114, 6236 (2001). 



11 



[28] C. Giardina and R. Livi, J. Stat. Phys. 91, 1027 (1998). 
[29] G. Rickayzen and J. G. Powles, J. Chem. Phys. 114, 4333 
(2001). 

[30] C.W. Gardiner and P. Zoller, Phys. Rev. A 58, 536 (1998). 
[31] C.W. Gardiner and M. J. Davis, cond-mat/0308044 (2003). 
[32] K. Binder, Z. Phys. B 43, 119 (1981). 

[33] P. Arnold and G. D. Moore, Phys. Rev. E 64, 0661 13 (2001). 
[34] L. Caiani, L. Casetti and M. Pettini, J. Phys. A 31, 3357 (1998). 
[35] L. Caiani, L. Casetti, C. Clementi, G. Pettini, M. Pettini and R. 



Gatto, Phys. Rev. E 57, 3886 (1998). 
[36] M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi and E. 

Vicari, Phys. Rev. B 63, 214503 (2001) 
[37] P. Arnold and G. Moore, Phys. Rev. Lett. 87, 120401 (2001). 
[38] V. A. Kashurnikov, N. V. Prokof 'ev, and B. V. Svistunov, Phys. 

Rev. Lett. 87, 120402 (2001). 
[39] G. Moore and P. Arnold, private communication (2003). 



